% This m-script estimates peer, contextual and individual effects
% Date: 2019-03-14
% Outputs:
% -- s1: row sum of reduced-form coefficients
% -- s2: column sum of reduced-form coefficients

function [BIND,BCLS,BCON,rBCON,rcm,rcmInd,rcmCls] = estSTAR(Yt,clsXt,indXt)
global varIndiv varClass dobwt LB UB meanVI meanVC meanY options inds
%% two-step estimation of social effects
for i = 1:2 % loop in the dummy 1('dobDistrn'<=6)
    for j = 1:2 % loop in small vs large classes
        %% reduced-form regression (pooled OLS using all classes)
        lb = LB(i,j); ub = UB(i,j);
        rcmInd{i,j} = zeros(lb,ub,length(varIndiv)); 
        rcmCls{i,j} = zeros(lb,length(varClass)); 
        for ii = 1:lb, % loop in individuals on l.h.s. of regression 
            rcmCls{i,j}(ii,:) = regress(Yt{i,j}(ii,:)', clsXt{i,j}');         
            for ji = 1:ub, % loop in classmates on r.h.s. of regression
                rcmInd{i,j}(ii,ji,:) = regress(Yt{i,j}(ii,:)', squeeze(indXt{i,j}(:,ji,:))');                 
            end;
        end;
        % calculate row sum and column sum of reduced form coefficients
        for kk =1:3,
            rcm(i,j,kk,:,:) = rcmInd{i,j}(1:15,1:15,kk);  
        end;        
        %% Find (a_k,b_k) using the implications of trace and column sums        
        mm{i,j} = zeros(length(varIndiv),1); 
        tt{i,j} = zeros(length(varIndiv),1); 
        for ik=1:length(varIndiv),
            M{i}(ik,j)  = mean(rcmInd{i,j}(:,:,ik)*ones(size(rcmInd{i,j},2),1));
            tt{i,j}(ik) = sum(diag(rcmInd{i,j}(1:lb,1:lb,ik))); 
        end;
        ab{i,j} = zeros(length(varIndiv)-1,2);
        for ik = 1:length(varIndiv)-1,
            ab{i,j}(ik,:) = inv([tt{i,j}(ik) tt{i,j}(end); M{i}(ik,j) M{i}(end,j)])*...
                                [lb;1];
        end;        
    end;
    A{i} = [ab{i,1}(:,1),ab{i,2}(:,1)];
    B{i} = [ab{i,1}(:,2),ab{i,2}(:,2)];
    %% Structural estimates using exclusion restrictions 
    % We maintain that 'days of absence' has no contextual effect; and that
    % -- peer effects vary with group size
    % -- direct & contextual effects do not var with group size
    [Mat,Vec] = constructLS(A{i},B{i},M{i});
    colInd = [1:5 7:8]; % column index for non-zero parameters
    fH = @(b) (Mat(:,colInd)*b-Vec)'*(Mat(:,[1:5 7:8])*b-Vec);
    index = nchoosek([1:14],length(colInd));
    ft = 1e2; fval = 1e2; b = zeros(7,1);
            for ri = 1:size(index,1), % ri,
                rows = index(ri,:);
                if rank(Mat(rows,colInd)) < length(colInd),
                    continue;
                else ...
                   est   = Mat(rows,colInd)\Vec(rows);
                   if (est(1)>.2)*(est(1)<.9)*(est(2)>.2)*(est(2)<.9)*...
                           (est(3)<0)*(est(4)>0)*(est(5)>0)==1; % *(est(6)>0)*(est(7)>0)==1,
                      [x,fval,exitFlag] = fminsearch(fH,est,options);
                   end;
                end;
                if fval < ft & exitFlag == 1,
                    ft = fval; b = x;
                end
            end;       
    bIND(:,i) = b; 
    %% structrual estimates for class variables
    bCLS(i,1) = mean(rcmCls{i,1})*(1-b(1));
    bCLS(i,2) = mean(rcmCls{i,2})*(1-b(2)); 
end; % bIND', bCLS, 
%% aggregate estimates 
BIND   = bIND*[dobwt;1-dobwt];  % BIND(1)=min(bIND(1,:)); BIND(2)=max(bIND(2,:));
BCLS   = bCLS'*[dobwt;1-dobwt];

%% estimates for structural intercepts
for i = 1:2, 
    for j = 1:2,
        ss(i,j) = 0; 
        for is = 1:UB(i,j)-LB(i,j)+1,
            for ii = 1:LB(i,j),              
                ss(i,j) = ss(i,j) + ...
                meanY{i,j}(is,ii) ... 
                 - sum(sum(squeeze(meanVI{i,j}(is,:,:).*rcmInd{i,j}(ii,:,:)),2))...
                    - meanVC{i,j}(is)*rcmCls{i,j}(ii);
            end;
        end;
        bCON(i,j) = (1-bIND(j,i))*ss(i,j)/(UB(i,j)-LB(i,j)+1)/LB(i,j);
    end;
end;
BCON = bCON'*[dobwt;1-dobwt];
rBCON = BCON./(1-BIND(1:2));
